Boolean Dynamics of Kauffman Models with a Scale-Free Network 
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We study the Boolean dynamics of the "quenched" Kauffman models with a directed scale- free 
network, comparing with that of the original directed random Kauffman networks and that of 
the directed exponential-fluctuation networks. We have numerically investigated the distributions 
of the state cycle lengths and its changes as the network size A*' and the average degree (fc) of 
nodes increase. In the relatively small network (A'^ ~ 150), the median, the mean value and the 
standard deviation grow exponentially with A'^ in the directed scale- free and the directed exponential- 
fluctuation networks with (fc) = 2, where the function forms of the distributions are given as an 
almost exponential. We have found that for the relatively large A'' ~ 10^ the growth of the median of 
the distribution over the attractor lengths asymptotically changes from algebraic type to exponential 
one as the average degree (fe) goes to (fc) = 2. The result supports an existence of the transition at 
{k)c = 2 derived in the annealed model. 
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I. INTRODUCTION 



The origin of life has attracted many scientists as one 
of the unsolved problems in science for a long time 1] . To 
answer the quest, the self-organization of matter^^l and 
the emergence of order Q have been regarded as the key 
ideas. Kauffman first introduced the so-called Kauffman 
model - a random Boolean network ^RBN) model, based 
upon the random network theory0|,(3,Q. This model has 
become a prototype for many authors to study complex 
systems such as metabolic stability and epigenesis, ge- 
netic regulatory networks, and transcriptional networks 
^ 2.1 A A ITH Il^ lIslla s well as general Boolean net- 
works [bi Ip-Tjli 111 ITsL ITgt ■ neural networks poj and 
spin glasses pit 123. l2a | . 

In the RBN, we assume that the total number of the 
elements (i.e., nodes) N and the link number of an ele- 
ment (i.e., the degree of a node) if in a directed random 
network are kept to be constants, respectively. Kauff- 
man found numerically that there is a phase transition 
from an ordered phase to a chaotic phase through the 
critical point (i.e., edge of chaos) sA K = = 2 for 
the equally probable model of the Boolean functions for 
and 1. Here in the chaotic phase the number of the 
state cycle attractors grows as exponential in N; in the 
ordered phase the growth is proportional to N; in the 
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critical point the growth is proportional to vA^. Later, 
this phenomenon has been analytically verified by study- 
ing the "annealed" Kauffman models 0, 0, ^^'^ 
also exact results have been analytically obtained for the 
special cases of if = TVf^llSilMEillllii and iC = 1 

Recently there has appeared a revival interest on the 
study of the critical point of the Kauffman models with 
Kc = 2 as well as = 1/p where p — 2p{l—p) andp de- 
notes a probability (0 < p < 1) given in the next section 
lilllllallllllllMlliliSlMllllil. Here, thepower- 
law distributions of the cycle length of attractor (33, the 
models with the finite numbers of 1 < if < 4 35, 3^.13^. 
the special cases of Kc = 2 networks such as the one- 
dimensional network and the Cayley treejs^, the scaling 
nature of the critical (ifc = 2) point and the chaotic 
{K > 2) phasefH, the superpolynomial growth of the 
number of the state cycle attractors ,40J, the stability of 
the system jij have been studied both numerically and 
analytically. These results suggest that the critical point 
of Kauffman model with if ^ = 2 is very special. 

Especially, BastoUa and Parisils^ |3^ Is^l elucidated 
that what is important on the critical line is the effective 
connectivity ife// (~ 1) and that the critical point of 
the model exhibits the nature of percolation of informa- 
tion where the character shares much with that of if = 1 
studied by Flyvbjerg and Kjaeri30||. To stu dy t his rather 
peculiar aspect further. Coppersmith et. al.^^H^I stud- 
ied numerically the "reversible" Kauffman model that 
is different from the original one in the sense that the 
former is dissipative while the latter is non-dissipative. 
They supplied a critical value of ifc = 1.62 for the non- 
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dissipative cases. 

On the other hand, at nearly the end of 1990's the 
scale-free networks have been discovered from study- 
ing the growth of the internet geometry and topology 
Qjt I451 11^ Isol Isij . After the discovery, scientists have 
become aware that many systems such as those which 
were originally studied by Kauffman as well as other var- 
ious systems such as internet topology, human sexual re- 
lationship, scientific collaboration, economical network, 
and so on. be long to the category of the scale-free net- 
works [il. I45I liaTlsol IsH . Therefore, it is very natural 
to apply the concepts of scale-free networks to the Kauff- 
man model and to ask whether or not there appears the 
difference between the RBN and the scale-free random 
Boolean network (SFRBN). 

Recently, as to this direction, the dynamics of the 
" qu enched" SFRBN has been intensively studied [s^, IH, 
m m m |5a |59l led, m m m. it has been shown 
that the critical point occurs when the average degree (fc) 
of the scale-free networks equals (fc)c = 2 for the equal 
probability (i.e. no bias) models of p = 0.5j_analytically 
studying the "annealed" SFRBN [El |^ |^, £QJ that 
are the scale-free analogs of the "annealed" RBN stud- 
ied by Derrida and Pomeaup^. Analytical results have 
been obtained for the annealed model with N —^ 00 un- 
der some assumptions such as mean field approximation 
and ergodicity. On the other hand, a lot of numerical re- 
sults have been calculated for finite-size quenched model. 
In the numerical simulation there are various ways for 
the sampling over, for example, (A) different networks, 
(B)different initial states, and (C)different Boolean func- 
tions assigned for each node. In particular, the mean 
value over the numerical distribution becomes insignifi- 
cance when the distribution is broad because the order 
of fluctuation exceeds over the mean value. Then the 
relative fluctuation diverges for system size going to in- 
finity, i.e. not self-averaging. Such a phenomena can be 
found at second-order phase transition and scale-free net- 
works with power-law correlation [s^ . In the situation 
the median value is statistically more robust against the 
artifacts due to undersampling of the set of the networks 
and the initial states than the mean value. 

In this paper we would like to mainly investigate a 
qualitative transition in the function form of the me- 
dian value fh over the distribution of attractor length as 
a function of the finite N, changing the average degree 
(k). Note that the transition does not always mean phase 
transition in the usual statistical mechanics sense, which 
is defined in the thermodynamic limit of A'^ — > 00. One of 
our interests is also finite-size effect of the transition phe- 
nomena around (fc) = 2 in the Boolean dynamics. How 
is the transition observed in the relatively small network 
corresponding to the realistic network size? 

To ease the reference, we first summarize our main re- 
sult: There is a transition of the function form m oc N" 
from a < 1 to a > 1 within a range of 1.4 < {k)c < 1-7 
in the "quenched" SFRBN without bias. The transi- 
tion becomes important when we consider some real- 



istic biological networks, as Kauffman first introduced, 
because the realistic data suggest that the number of 
cell types in organism is crudely proportional to the lin- 
ear or square-root function of the DNA content per cell, 
(X N°'{a < I). We adapted a completely synchronous 
updating rule in the Boolean dynamics because the at- 
tractor lengths are well-defined only for the synchronous 
updating rule. However, it is noting that the synchronism 
idealization is not always true for biological systems as 
genetic regulatory networks and asynchronous updating 
rules are more plausible for biological phenomena |64j . 

Moreover, the function form fh{N) asymptotically 
changes from the algebraic type m{N) oc A^" to the ex- 
ponential one as the average degree (fc) goes to (fc) = 2. 
This result is consistent with the previous belief that the 
transition occurs at the critical value of the degree of 
nodes of ATc = 2 in the "annealed" RBN and SFRBN 
si 

The organization of the paper is the following. In 
Sec. II, we present the Kauffman models that we study. In 
Sec. Ill, to apply various networks to the Kauffman mod- 
els we give how to generate such networks. In Sec. IV, we 
study the cycle distributions of the " quenched" Kauffman 
models for the various network systems. In Sec.V, con- 
clusions will be made. We mainly focused on the change 
of the functional form of the median around the critical 
value (fc) < 2 in the main text. However, the other quan- 
tities such as Derrida plots and frozen nodes density are 
also used in order to investigate the Boolean dynamics. 
We briefly give analytical result for (fc) = 1 in appendix 
A and Derrida plots in appendix B. 

II. RANDOM BOOLEAN NETWORK (RBN) 

The RBN requires us to assume that the total num- 
ber of nodes (vertices) A^ and the degree (the number 
of links) of the i-th node ki are fixed in the problem, 
where all ki = K . Since there are K inputs to each node, 
2^ Boolean functions can be defined on each node; the 
number 2^ certainly becomes very large as K becomes 
a large number. Then, we assume that Boolean func- 
tions are randomly chosen on each node from the 2^ 
possibilities. Locally this can be given by 

a,{t + l) = B,{<j,,{t),a,^{t),--- ,<T,,^(t)) (1) 

for i = 1, . . . , A^, where e Z2 = {0, 1} is the binary 
state and Bi e Z2 is a Boolean function at the ith node, 
randomly chosen from 2^ ' Boolean functions, where the 
probability to take 1 (or 0) is assumed to be p (or 1 — p) 
(TABLE 1). In this paper we used a case oi p — 1/2 
for the numerical calculations. But the treatment for 
other cases of p 7^ 1/2 is straightforward[65j . If we 
fix the set of the randomly chosen Boolean functions 
{Bi,i = 1, • • • , A^} in the course of time development, 
then this model is called "quenched' model^]. On 
the other hand, if we change each time the set of the 
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randomly chosen Boolean functions in the course of time 
development, then this model is called " annealed' model 
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TABLE 1. The relationship between the Boolean 
functions Bi and inputs, {<7ij , cr^j , • • • ,<Ji^ }. Since there 
are ki Ci's, each of which has or 1, there are 2'^' ways of 
inputs. These provide 2''^ai's of outputs, each of which 
is or 1, randomly chosen with a probability of p or 
1 — p. Hence, there are possibly 2^ ' Boolean functions 
at each node. 



we obtain the critical value pc = (1± -^1 — 2/K)/2. As 
is described in the introduction, the analytical methods 
for the systems with special values of the degree of nodes 
'27,'25','M,|2l|2i|2a|33,|3l|43 have been studied in 
details, already. 



III. VARIOUS NETWORK MODELS 

To apply the RBN to that with a given network, 
we have to specify what kind of network model we 
take. Since there are so many types of networks |4^. we 
limit ourselves only to consider the directed random net- 
works that were first considered by KauffmanQ, the di- 
rected scale-free networks, and the directed exponential- 
fluctuation networks in this paper. But the generaliza- 
tion to other networks is straightforward. In this section, 
we will give how to generate such directed networks ex- 
cept the directed random networks, since the generation 
of those is very well known 0, 0, . 



We then study the dynamics of the "quenched" RBN. 
Since there are N nodes in the random network, there 
are 2^ states jV') = jci, a2, ■ ■ ■ , cat) in the system such 
as iVo) = |0,0,0,...,0),...,|V'2«-i) = 
These states form all vertices of a hypercube in N dimen- 
sions. First, we start to define the initial state IV'^^'') out 
of the 2^ states (for example, say = |0, 1, 1, . . . , 0)). 

Second, according to the Boolean functions defined on 
the network, the initial state evolves to another state 
lip^"^^) in the 2^ states. Hence, the time development can 
be simply denoted as 

= B.I^W), (2) 

where each entry in can be developed by Eq.(l) and 
Bt means a Boolean function chosen at time t. Therefore, 
if we collect all the 2^ states as a 2 ^-dimensional vector 
|<If(t)) = (|0, 0, 0, . . . , 0), . . . , |1, 1, 1, . . . , 1))^, then we can 
get a matrix representation of Eq.(2) as 

= Bt\¥'^). (3) 

Here Bt is a 2^ x 2^ matrix that all components are only 
or 1, and means a state vector whose compo- 

nents are degenerate such that the mapping of Eq.(3) is 
not a one-to-one but a many-to-one correspondence. 

Considering this time developing equation, we find the 
cyclic structure of the states such as the length of the 
state cycle, the transient time and the basin size, and 
so ony, Qj It is obviously difficult to do this pro- 
cedure by an analytical method in general. Therefore, 
we must do it numerically,!^ as well as analytically if 
possible. However, there are some important analytical 
results. The analytical investigations on the "annealed" 
models 0, 16, 34] showed the existence of a phase 
transition at the critical value of Kc = l/[2p(l — p)] 
[T^ ITEl ll^ l3^ , and if we solve it conversely for p then 



A. Directed scale-free networks with the integer 
average degree of nodes 

Let us first consider the directed scale-free networks. 
In this case, we adopt a little modified version of the so- 
called Barabasi- Albert model 45], since we have to treat 
the directed scale-free networks with fractional numbers 
of the average degree of a node (i.e., vertex), (fc), such as 
1 < (fc) < 2. 

Denote by Vi the i-th node to which we want to put 
links and denote by Vj one of the surrounding nodes. 
Then, the input from the j-th node to the i-th node is 
described by the in-going arrow as Vi <— Vj, while the 
output from the i-th node to the j-th node is described 
by the out-going arrow as Vi Vj. Denote by /c™ (= 
kj^i) the input degree of the i-th node. Denote by fc""* 
(= fcjv-i) the output degree of the i-th node. We note 
that from simple consideration, when one input link to 
the i-th node is put between the i-th node and the j-th 
node, it becomes one output link for the j-th node at the 
same time. Therefore, if one input link is increased, so is 
one output link simultaneously. 

Let us consider the case that the average degree of 
nodes, (fc), is an even integer, such as (fc) = 2n. (Al) 
We initially start with mo (> n) nodes for seeds of the 
system. We assume that both input and output links 
are simultaneously linked between all of the initial seed 
nodes. Therefore, the total link numbers for the input 
and the output are n{n — l)/2 = Lq, respectively. (A2) 
Next, every time when we add one node to the system, 
m ~ n {< mo) new links are randomly chosen in the 
previously existing network, according to the preferential 
attachment probability for the output network. 

Lout 
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(A3) Similarly we redo the same procedure for the input 
network, according to the preferential attachment prob- 
ability for the input network, 



n,(fcf) 



(5) 



We continue the above procedures until the system size 
N is achieved. Therefore, after t steps, we obtain the 
total number of nodes as N{t) = toq + t and the total 
numbers of links for the input U"{t) and output 
as = L°"*(i) = Lq + 2nt, respectively. Hence, 

by this we can obtain for the directed scale-free network 
(fc™) EE 2L'''{t)/N{t) = 2n, EE 2i°"*(t)/iV(t) = 2n 

as i ^ oo. 

When (fc) — 2n + 1, after the procedure (A3) we add 
one more procedure: (A4) we redo the procedure (A2) or 
(A3) with equal probability. 



B. Directed scale- free networks with tlie fractional 
average degree of nodes 

Let us next consider the directed scale-free networks 
with a fractional average number of degree. Suppose that 
(fc) is fractional such that (k) = [(fc)] -I- {{k}) where [{k}] 
means the integer part (say, n) and ((fc)) the fractional 
part (say, /) such that {k) = n + /, where < / < 1. 

In this case, (Bl) we first follow the same procedure 
(Al) in the previous subsection A. (B2) Second, we add 
one node to the system at each step of time. Every time 
when a new node is added, we have to define the node to 
give input or output. For this, we randomly choose in- 
put or output with equal probability of 0.5. (B3) Third, 
if the chosen case is input (output) for the node, then 
we follow the procedure (A2) [(A3)] in the previous sub- 
section A. Then, we place the input (output) links with 
equal probability 1 among the chosen links. (B4) Fourth, 
if the not-chosen case is output (input) for the node, then 
we follow the procedure (A3) [(A2)] in the previous sub- 
section A. Then, we place the output (input) links with 
equal probability //n (0 < / < 1) among the chosen 
links. (B5) Fifth, go back to (B2) and redo the same 
procedures until the system size N is achieved. 

Then, after t steps, we obtain the total number of 
nodes N{t) = n + t, the total number of input links 
U'^{t) = io + ("■ + /)^, and the total number of out- 
put links L°^*{t) — Lo + {n + f)t, respectively. Hence, we 
can obtain the average input and output degrees of nodes 
(fc'") = U''{t)/N{t) = n+f and (fc°"*) = L'"'\t)/N{t) = 
n + f as i ^ oo , respectively. 

Using the above method, we can construct a scale-free 
network with a fractional average degree of nodes 1 < 
(k) < 2. For example, consider the case of (fc) — 1.4. In 
this case, we just take n = 1 and / = 0.4. 






J-r 




FIG. 1: (Color online) Typical examples of directed net- 
works are shown for the size of TV = 64. (a),(b) Random 
network with K — 2; (c),(d) Exponential-fluctuation network 
with (fc) = 2; (e),(f) Scale-free network with (k) = 2. We 
show same network by two kinds of representation. For (a), 
(c) and (e) the nodes are located on the circumference with 
equal distance. For (b), (d) and (f) the nodes are randomly 
distributed in the square. Each node is represented as a bold 
point with size in proportion to the number of the input links. 
We represent input (output)-side of the links with deep(faint) 
color such that the direction of a link is denoted by the color 
gradation from deep color(output) to faint color(input). 



C. Directed exponential-fluctuation networks 

Let us consider the directed exponential-fluctuation 
networks. We can follow the same procedure for both 
the cases of the integer and fractional average degrees of 
nodes, replacing the probabilities of Eqs.(4) and (5) by 



H,(fcr)=n,(A:r*) = 



1 



mo + 1 



(6) 



The exponential-like distributions are often observed in 
some real- world networks • The generalization to 
other networks can be straightforward 49j . 

Fig.l shows typical examples of the directed networks 
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where the total number of nodes A'' = 64 is used: (a),(b) 
the directed random network with the degree of nodes, 
K — 2; (c),(d) the directed scale-free network with the 
average degree of nodes, (fc) — 2; and (e),(f) the directed 
exponential-fluctuation networks with the average degree 
of nodes, (fc) = 2, respectively. 



IV. CYCLE DISTRIBUTIONS 

Now we are going to apply the above-mentioned 
"quenched" Kauffman's Boolean dynamics to the di- 
rected random networks, the directed scale- free networks, 
and the directed exponential-fluctuation networks. The 
first one provides the famous "quenched" RBN, the sec- 
ond one the "quenched" SFBN model, and the third one 
the "quenched" exponential- fluctuation random boolean 
network(EFRBN) model 1^. 

Before going to present the numerical results, let us 
explain the calculation method for obtaining the lengths 
of the state cycles as follows: (i) Realizations: We first 
prepare 1,000 sample networks with the size of N and 
investigate them in order, (ii) Initial conditions: Second, 
we randomly choose an initial state out of the 2^ 

states, (iii) The "quenched" Boolean dynamics: Third, 
we calculate the "Boolean dynamics" [given by Eq.(2)] 
repeatedly from 10 times to 100, 000 times, where we fix 
a special set of the Boolean functions on the network. 
This means that we treat the "quenched" model, (iv) 
State cycles: Fourth, we investigate whether there ex- 
ist the states that belong to the equivalent state in the 
data. Since the "Boolean dynamics" is deterministic, if 
the state |vl/(*i)^ at time t — ti is the same as the state 
|^(ti+s)^ at t = ti+s, then the state at t = ti + 1 

becomes the same state as the state |v[/(*i+^+i)) at time 
t = ti + s + 1. Therefore, the length of the state cycles 
satisfying this condition is s -f 1. Such calculations are 
performed for all network samples. 
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FIG. 2: (Color online) Histograms h of the lengths £c of state 
cycles in various types of the directed networks. The network 
size is TV = 40. Each histogram is generated by 10^ differ- 
ent sets of the Boolean functions and five different network 
structures. The maximum iteration number of the Boolean 
dynamics is 10^ until the convergence to the cycle is realized. 
(a)-(b) RBN with K = 2,4- (c)-(d) EFRBN with (fc> = 2,4; 
(e)-(f) SFRBN with (fc> = 2,4. Here in (a)-(b) 2^"^ Boolean 
functions of K variables are used equiprobably, while in (c)-(f ) 

2^ ' Boolean functions of ki variables are used equiprobably 
at the i-th node. The lines are the exponential fittings to the 
numerical data. The mesh sizes of the histograms are 1 for 
A' = 2 and (k) = 2, 10^ for K = 4 and (k) = 4. 



B. The median fh, the mean value (Ic), and the 
standard deviation a of the distributions of the state 
cycle lengths 



A. Histograms of the lengths £c of the state cycles 

Fig. 2 shows the histograms of the lengths £c of the state 
cycles in the original "quenched" RBN with K = 2,4 and 
in the "quenched" SFRBN and EFRBN where the aver- 
age degree of nodes (fc) ~ 2, 4. We can see that the func- 
tional forms of the distributions seem exponential type. 
The tail of the distribution depends on the fluctuation 
property of the degree of the nodes. We found that in 
comparing among the three types of the network struc- 
tures the maximum length of the state cycle becomes 
longer as the fluctuation in the degrees of nodes becomes 
larger. We try to investigate the more accurate functional 
form of the distribution and the transition depending on 
the degree of nodes K and (fc). 



Fig. 3 shows the median fh, the mean value {£c), and 
the standard deviation a of the distributions of the state 
cycle lengths with respect to the total number 150) 
of nodes for the various directed networks, respectively. 
Fig. 3(a) shows that in the RBN with K — 2 the median 
grows as \/N in relatively small N region and it grows 
as proportional to N in the large N. Whether or not the 
behavior of -v/lv is valid is very delicate since we have to 
always stem this from the numerical data of the finite size 
systems. The more details will be shown elsewhere [gs). 

On the other hand, in the RBN with K = 4 the median 
grows exponentially with as m ^ e"^. We have ob- 
served that as K increases, the growth type of the median 
exhibits a transition from algebraic type to exponential 
one in the RBN. In the EFRBN and the SFRBN, the 
median grows exponentially with N for (fc) = 2 and 4, 
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FIG. 3: (Color online) Semilog-plots of the the median fh, 
the mean value {£c}, and the standard deviation a of 5000 
samples of the state cycles with respect to the total number 
of nodes for the various directed networks, respectively. 
(a)-(c) are shown for the RBN with K = 2,4. (e)-(f) are 
shown for the EFRBN with (k) = 1,2,4. And (g)-(i) are 
shown for the SFRBN with (fc) = 1,2,4. The lines are best 
exponential fittings, m ~ e"^, {£c) ~ e*^, a ~ e'^^ , to the 
numerical data except for the RBN with K = 2. The best 
fitted exponents are a = 0.02, b = 0.037, c = 0.044 for the 
EFRBN with (k) = 2; a = 0.0168, b = 0.022, c = 0.025 for 
the SFRBN with {k) = 2; a = 0.163,6 = 0.193, c = 0.194 for 
the RBN with if = 4; a = 0.170, b = 0.179, c = 0.194 for the 
EFRBN with {k) = 4; a = 0.123, 6 = 0.131, c = 0.136 for the 
SFRBN with (k) = 4. 
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FIG. 4; (Color online) The relationship between the mean 
value (Ic) and the standard deviation a of the distributions of 
the state cycle lengths with respect to the total number A'' of 
nodes for the RBN with K = 2,4 and the SFRBN {k) = 2, 4, 
respectively. The lines are fitting to the numerical data as 
a ~ (ic)'^. The best fitted exponents are (3 = 1.61, 1.02 for 
the RBN with if = 2,4, and f3 = 1.19, 1.02 for the SFRBN 
with (k) = 2, 4, respectively. The resuh for the EFRBN is 
almost the same as that for the SFRBN. 



The relationship between the mean value (ic) 
and the standard deviation a 



Fig. 4 shows the relationship between the mean value 
(ic) and the standard deviation a of the distributions of 
the lengths of the state cycles with respect to the total 
number N of nodes for the RBN and the SFRBN, respec- 
tively. We have found that the relationship is fitted by 
the following expression: a ~ (4)^. The best fitted ex- 
ponent /9 is nearly equal to unity for K = A and (k) = 4. 
As a matter of fact, we can guess that in the continuum 
limit of the relatively large TV, the distribution P{£c) ap- 
proaches the exponential form P{ic) = Q;exp(— a£c), as 
K and (fc) increases. In the exponential distribution, the 
mean value (ic) and the standard deviation a can be 
represented in terms of the single parameter a such as 
(4) =<7^ 1/a. 



and such a transition can not be observed in the range 
between (k) = 2 and (k) — 4. Therefore, we conjec- 
ture that based on Fig. 3, the transition takes place in 
the range between (fc) = 1 and (k) = 2 in the RBN and 
the EFRBN. Note that N dependence for (fc) = 1 can be 
analytically derived. (See appendix A.) 

Furthermore, we can see the mean value {£c) and the 
standard deviation a grow exponentially with N such 
as {£c) ~ e''^ , a ~ e'^^, in all cases. It is difficult to 
distinguish clearly the dynamical properties between the 
SFRBN and the EFRBN when the network size N is as 
small as 150. 



D. Plots of the median value m 

As shown in Fig. 4 it is clear that the standard devi- 
ation (T is larger than the mean value {£c) ioi K — 2 
and (k) — 2. In the cases, the numerical mean values are 
sometimes without credibility due to the lack of the sam- 
ples with large cycle lengths. Such an imperfection of the 
self-averaging is well-observed for the distributions with 
a long-tail. Therefore, we investigate the relatively stable 
median as representative values to characterize the dis- 
tributions of the attractor lengths, instead of the mean, 
for the large networks. 

Fig. 5 shows the median value m of the distribution of 
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FIG. 5; (Color online) Plots of the median value m of 1000 
samples of the lengths of the state cycles with respect to the 
total number A*" of nodes for the directed scale-free networks 
with (fc) =1.3,1.4,1.5,1.6 and 1.7. The results for iV up to 1024 
are shown. We used a cut-off attractor cycle length lent ~ 10^ 
due to the numerical difficulty to detect the large cycle length. 
For the 10-20 network samples, we did not detect the attractor 
cycle lengths. However, it seriously does not influence on the 
median of the distribution. 



state cycle lengths with respect to the total number N 
of nodes for the relatively large SFRBN {N - 10^). We 
found that there is a transition in a range of 1.4 < {k)c < 
1.7, dividing the N"{a < 1) growth and N°'{a > 1) 
growth in polynomial growth in the ordered phase. Note 
that some biological data suggest slow growth as N" {a < 
1) in the relatively small network size {N ~ 10'^). 

It is well-known that the scale- free topologies are ubiq- 
uitously in nature and the degree exponents lie in be- 
tween 2 and 3. Furthermore, one of the important prop- 
erty in the scale-free topology is the existence of the 
highly connected hub as seen in the yeast synthetic net- 
work and so on [b^ . The realistic systems do not contain 
enough nodes to closely approximate the true transition, 
therefore, the finite-size behavior is important. In fact 
some biological realistic data suggest that the number of 
cell types in organism is crudely proportional to the lin- 
ear or square-root function of the DNA content per cell, 
i.e. rh cx N°'{a < 1) 3], for the finite-size TV ^ lO'^. 

Finally, we investigate a transition between the poly- 
nomial growth in the ordered phase and the exponential 
one in the chaotic phase. Figure 6 shows the semi-log and 
log-log plots of the data near (fc)2 '--^ 2 in Fig. 5. It is sug- 
gested that the function form m{N) has an polynomial 
form for (fc) < 2 and approaches the exponential form as 
{k) ^ 2 in the quenched SFRBN. The results are con- 
sistent with the critical value Kc — 2 given for annealed 
SFRBN in other reports [HililH- In this section 
we focused on the scaling of attractor length near criti- 
cal value (fc) ^ 2. The other quantities are also useful 
measure for the transition of the Boolean dynamics. (See 
appendix B.) 




N 

FIG. 6: (Color online) (a)Semilog-plots and (b) log- log-plots 
of the median m with respect to the total number of nodes 
for some values (k) near (fe) = 2 for the scale-free network 
in Fig. 5. Apparently, the case with (fc) — 2 exponentially 
increases with respect to A^. The straight line in (a) shows 
the exponential growth for a guide for eyes, and the straight 
lines in (b) show the polynomial growth (m oc N") with the 
slopes a — 2.5, 1.8, 0.98 from top to bottom. 

V. CONCLUSIONS 

In conclusion, we have studied the Boolean dynamics of 
the " quenched" Kauffman model with directed scale- free 
networks(SFRBN), comparing with ones of the directed 
random networks (RBN) and the directed exponential- 
fluctuation networks (EFRBN). We have numerically cal- 
culated the distributions of the lengths of the cycles and 
its changes as the network size N and the average degree 
of the node increase. We have found that the median, 
the mean value and the standard deviation grows expo- 
nentially with N in the EFRBN and the SFRBN with 
(k) — 2,4, where the function forms of the distributions 
are almost exponential. 

/^From our results we conclude that a transition occurs 
near (k) ^ 2 in the SFRBN. The result is consistent with 
that in the "annealed" SFRBN. 

In this paper we dealt with the directed graphs that 
correspond to the asymmetric adjacency matrices in the 
network theory. However, the quite different distribu- 
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tions of the state cycle lengths are observed in the undi- 
rected random networks. Therefore, it is also interesting 
to study the properties of the distributions of the state 
cycle lengths in the undirected scale- free networks as well. 
The details will be given elsewhere [6^ . 

We have numerically investigated the finite-size net- 
works {N - 140 - 1024). In fact, it is very difficult 
to study the indefinitely large size of systems. This 
seems the disadvantage of our approach of the finite 
size systems. However, practically speaking, there are 
a lot of scale-free-type biological networks with size N ^ 
O(IO^) — 0(10'^), such as the metabolic reaction networks 
of the bacterium E.coli and of the Yeast protein interac- 
tion networks, etc0) EBl- Thus, we beheve that our 
results stemmed from the finite size systems might play 
an important role to study such real network systems in 
Nature. As we mentioned in the introduction, the asyn- 
chronous updating rules are more plausible for biological 
phenomena. It is interesting to expand the investigation 
to the asynchronous version. 



to the Lyapunov exponent in the continuous dynamics, 
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FIG. 7: (Color online) Mean value (£c) of attractor lengths 
in RBN and SFRBN with (fe) = 1. The analytical result 
{£c) = -p/(l-p) log(l-p) ~ 1.44 for p= 1/2 in the annealed 
RBN is also overplotted. 



Acknowledgments 

We would like to thank Dr. Jun Hidaka for collecting 
many relevant papers on the Kauffman model and the 
related topics. One of us (K. I.) would like to thank 
Kazuko Iguchi for her continuous financial support and 
encouragement . 

APPENDIX A: A CASE WITH (k) = 1 

In this appendix, we compare the numerical results in 
SFRBN with (k) = 1 with analytically result in RBN 
with K = 1. 

In Fig. 7, we show the finite-size effect of the mean 
value {£c}. The analytical result (4) ~ 1.44 in RBN 
with K — 1 have been derived based on the probability 
for distribution of information conserving loops by Fly- 
vbjerg and Kjaer |30j . It is found that the numerical 
results converge a value 1.44 as the system size increases 
in both SFRBN and exponential-fiuctuation network. 

As mentioned in Sect. IV, we constructed the distribu- 
tion of the attracter length by different networks with 
(k) = 1, i.e. case (A) in introduction. Then the states 
converge into a point attractor with period unity, £c = 1, 
through the transient states in a lot of the network sam- 
ples. We can partially see the result in the distribution 
in 121 Accordingly the median fh over the distribution is 
always unity independent of syatem size N in our case. 

APPENDIX B: DERRIDA PLOTS 

In this appendix, we give the numerical results of Der- 
rida plot in the SFRBN. The Derrida plot, analogous 



measures the divergence of trajectory based on the nor- 
malized Hamming distance between two distinct states 
a}{t),afit),i = l,2,..,N, 



Hit) = -Y.\al{t)~af{t)\. (Bl) 



For a sample of the random initial states pairs, the av- 
erage H{t + 1) is plotted against H{t), and the plot is 
repeated for increasing H{t). A curve above the main di- 
agonal H{t + 1) = H{t) indicates a divergent trajectory 
and chaotic. The area below the diagonal H{t+1) = H{t) 
is called the ordered region, because the trajectories con- 
verge in the state space. A curve tangential to the main 
diagonal H{t + 1) = H{t) indicates a critical dynamics. 

Figure 8 shows the Derrida plots corresponding to the 
SFRBN with (k) = 1, 2, 4. It is apparent that the SFRBN 
with large (k) exhibits more chaotic behavior than that 
with small (k). It is important to note that the SFRBN 
is more ordered than the RBN compared with the cases 
with K = (fc). The Derrida coefficient is defined by 
Dc = logs, where s is the slope of the Derrida plot 
(curve) at the origin. If s = 1 then Dc = 0, giving 
the critical evolution. The slope of the Derrida plot for 
the SFRBN with (fc) = 1 is approximately unity. The re- 
sult is consistent with the occurence of the transition at 
(fc) = 2 in the function form of fh{N) observed in Subsec. 
IV.D. 

The other quantities such as, the density of frozen 
nodes and the robustness against perturbation, are also 
important indicators for the dynamical behavior in the 
RBN and the SFRBN. 
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Hamming Distance iH(t) 

FIG. 8: (Color online) Derrida plots of the SFRBN with 
(k) = 1, (fc) — 2 and (k) = 4. The analytical curves for the 
RBN with K = 1, K = 2 and A' = 4 are also overplotted.QA 
line H{t + 1) = H{t) is the dividing line between order and 
chaos. It is clear that K — 2 lies directly on this line, the 
system size is A'^ = 1024 and the number of the initial states 
for averaging is 2000. 
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